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Statement  of  Work 

Random  variables  arising  in  communications,  radar,  photodetection,  and  sig¬ 
nal  processing  often  have  moment  generating  functions,  characteristic  functions, 
or  probability  generating  functions  for  which  formal  mathematical  expressions  are 
readily  derived  theoretically.  Determining  the  probability  distributions  of 
these  variables  from  generating  functions  analytically  Is  seldom  possible.  The 
cumulative  distribution  of  such  a  variable,  important  because  it  provides  false- 
alarm,  detection,  and  error  probabilities,  can  be  expressed  as  a  contour  integral 
in  the  complex  plane  whose  integrand  involves  the  generating  function.  It  has 
been  found  that  the  distribution  can  be  accurately  computed  numerically  if  the 
contour  is  taken  through  a  saddlepoint  of  the  entire  integrand.  We  propose  a 
thorough  study  of  the  potentialities  of  this  method  of  saddlepoint  integration 
for  distributions  of  both  continuous  and  discrete  random  variables  of  importance 
in  the  above-mentioned  fields.  When  the  random  variable  results  from  a  quadratic 
functional  of  a  Gaussian  random  process,  the  generating  function  involves  the 
Fredholm  determinant  of  the  autocovariance  function  of  the  process.  Its  calcula¬ 
tion  for  points  on  the  contour  of  integration  by  integrating  vector  Riccati  equa¬ 
tions  will  be  investigated  for  processes  that  can  be  modeled  as  the  output  of  a 
linear  system  driven  by  white  noise.  A  coherent  signal  component  of  the  input 
process  can  also  be  handled  by  this  analysis.  Acceleration  of  the  convergence  of 
the  numerical  integration  by  optimum  choice  of  the  contour  will  be  examined. 

Particular  applications  treated  in  the  proposal  are  to  the  detectability  of 
fading  radar  echoes,  the  distribution  of  the  average  power  of  a  Gaussian  random 
process,  the  performance  of  the  optimum  and  threshold  detectors  of  a  narrow-band 
Gaussian  stochastic  signal  in  white  noise,  the  distribution  of  the  number  of 
photoelectrons  counted  when  a  mixture  of  coherent  and  incoherent  light  falls  on 
an  emissive  surface,  and  the  distribution  of  the  number  of  output  electrons  from 
an  avalanche  photodiode.  Still  other  applications  will  be  considered  if  time  and 
resources  permit. 
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1 .  Introduction 


Let  x  be  a  continuous  random  variable  with  probability  density  function 
(p.d.f.)  p(x)  and  moment  generating  function  (m.g.f.)  " 


h(u)  =  E(e_ux) 


•/' 


p(x)  e  ^dx. 


(1.1) 


Then  the  p.d.f.  p(x)  can  be  recovered  by  the  inverse  Laplace  transform 

.c+i« 


p(x) 


■L 


nv 

h(u)  e  du/2ni, 


(1.2) 


in  which  the  contour  of  integration  runs  parallel  to  the  imaginary  axis  and  lies 
in  a  strip  containing  the  origin,  but  no  singularities  of  h(u).  In  many  problems 
in  communications  and  signal  processing  the  m.g.f.  h(u)  is  readily  determined 
theoretically,  but  calculating  the  p.d.f.  p(x)  from  it  analytically  may  be  diffi¬ 
cult  or  impossible.  It  is  usually  the  cumulative  distribution  (c.d.) 


■f 


q  (x)  -  /  p(x')  dx' 

or  the  complementary  cumulative  distribution  (c.c.d.) 


(1.3) 


q  (x)  *  1  -  q  (x) 


p(x')  dx' 


(1.4) 


that  is  of  Interest,  for  these  determine  quantities  such  as  the  false-alarm  and 
detection  probabilities,  or  the  average  error  probabilities  in  a  communication 
system.  The  necessity  of  integrating  p(x)  enhances  the  complexity  of  the  prob¬ 
lem.  We  call  q”(x)  and  q+(x)  the  "tail  probabilities"  and  concentrate  on  calcu¬ 
lating  the  former  for  x  <  E(x)  and  the  latter  for  x  >  E(x),  where  E(x)  stands  for 
the  expected  value  of  the  random  variable  x.  They  are  given  in  terms  of  the 
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Z'.i  t  c , 


t  1 


; !  n 


m.g.f.  h(u)  by  the  contour  integrals 

c  >  0,  (1.5) 

c  <  0.  (1.6) 

The  principal  purpose  of  this  research  is  to  investigate  methods  of  compu¬ 
ting  the  tail  probabilities  q“(x)  and  q+(x)  by  numerical  integration  of  (1.5)  and 
(1.6)  along  a  suitable  contour.  Ideally  one  would  like  to  integrate  along  that 
contour  on  which  the  magnitude  of  the  integrand  decreases  most  rapidly  from  its 
value  at  the  point  Im  u  *  0  as  the  point  u  moves  away  from  the  real  axis  in  the 
u-plane,  for  then  the  number  of  steps  required  would  be  minimum  [1].  This  con¬ 
tour  is  known  as  the  path  of  steepest  descent,  and  along  it  the  integrand  is 
real,  or  equivalently  the  imaginary  part  Im  #(u)  of  its  "phase" 

♦(u)  “  In  h(u)  +  ux  -  In  (i  u)  (1.7) 

vanishes.  (Here  (±  u)  =  -u  when  (1.6)  is  being  Integrated  and  (±  u)  =  u  when 
(1.5)  is  being  integrated.)  This  contour  crosses  the  real  u-axis  at  a  saddle- 
point  u  -  Uq  of  the  integrand,  which  is  determined  by  the  equation 

*'(u)  -  +  x  -  u  1  -  0,  u  *  u0  (1.8) 

(Primes  stand  for  derivatives.)  One  such  saddlepoint,  uq  *  Uq-,  lies  between  the 
rightmost  singularity  of  h(u)  in  Re  u  <  0  and  the  point  u  *  0;  another,  uq  ■  uq+, 
lies  to  right  of  the  origin,  but  to  the  left  of  any  singularities  on  Re  u  >  0. 
When  (1.5)  is  being  integrated,  the  contour  of  integration  passes  through  uq+; 
when  (1.6)  is  being  integrated,  it  passes  through  Uq~.  When  the  second 
derivative  *"(u)  of  the  phase  is  easily  calculated,  (1.8)  can  be  expeditiously 
solved  by  Newton's  method.  Otherwise  one  can  use  the  secant  method,  which 
estimates  the  second  derivative  by  computing  *'(u)  at  two  close  points.  An 
alternative  method  for  determining  the  saddlepoint  is  to  solve  the  equation 


/cn» 
u"1 

■i“ 

/c+i® 


h(u)  e11*  du/2iri, 


h(u)  e**  du/2  id. 
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(1.9) 


Im  #(Uq  +  ie)  =  0 

by  Newton's  method  or  the  secant  method  for  a  sufficiently  small  value  of  e.  It 
is  unnecessary  to  determine  the  saddlepoint  with  high  precision. 

To  trace  the  path  of  steepest  descent  exactly  would,  however,  much  protract 
the  numerical  integration.  When  the  singularities  of  the  integrand  all  lie  on  or 
near  the  real  u-axis,  the  path  of  steepest  descent  often  has  approximately  a 
■*  parabolic  form,  and  the  contour  integrals  can  be  efficiently  evaluated  by  taking 
them  instead  along  the  osculatory  parabola,  that  is,  along  a  parabola  passing 
through  the  saddlepoint  Uq,  lying  symmetrically  about  the  real  u-axis,  and  having 
at  the  saddlepoint  the  same  curvature  <  as  the  path  of  steepest  descent,  which  is 
given  by  [1] 

x  =  j  *"'(uo)/*,,(u0).  (1.10) 

The  c.d.  and  the  c.c.d.  are  then  given  by  integrals  of  the  form 

q*(x)  =  /  Rl[(h(u)  e™  (+  u)  1  (1  -  ixy)]dy/it, 

J0  (1.11) 

2 

u  -  uQ  +  V2  *y  +  iy» 

with  Uq  =  uQ+  for  q— (x)  and  Uq  =  Uq~  for  q+(x) ,  and  for  reasons  given  in  [2]  this 

has  been  evaluated  by  the  trapezoidal  rule.  One  takes  an  initial  step  size  of 

-1  /2 

the  order  of  [$"(uq)]  ,  halving  the  step  size  and  repeating  the  numerical 

integration  until  the  value  of  the  probability  stabilizes  within  the  number  of 
significant  figures  desired. 

If,  on  the  other  hand,  the  integrand  also  possesses  singularities  above  and 
below  the  real  u-axi3,  and  the  path  of  steepest  descent  consequently  possesses 
many  branches,  it  may  not  be  feasible  to  integrate  along  other  than  a  straight 
line,  possibly  vertical,  possibly  inclined  at  a  suitable  angle  to  the  vertical. 

In  what  follows  we  shall  describe  -our  progress  in  this  Investigation,  treat¬ 
ing  in  Sec.  2  the  evaluation  of  radar  detection  probabilities  and  in  Sec.  3  the 
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calculation  of  the  cumulative  distribution  of  the  average  power  la  a  Gaussian 
random  process.  Similar  methods  can  be  used  for  computing  the  cumulative  distri¬ 
bution  of  a  positive  integer-valued  random  variable,  and  these  are  described  in 
Sec.  4  and  applied  to  distributions  of  the  number  of  electrons  at  the  output  of  a 
photomultiplier. 
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2.  Quadratic  Detection  of  Radar  Signals. 

(a)  Signals  of  Fixed  or  Random  Strength 

The  input  to  a  radar  receiver  is  passed  through  a  filter  matched  to  the 
shape  of  the  expected  echo  signal,  and  the  output  of  the  filter  is  applied  to  a 
quadratic  rectifier.  The  noise  in  the  input  is  white  and  Gaussian;  echoes  from 
a  target  may  also  be  present.  The  sampled  outputs  of  the  rectifier  in  a 
particular  range  bin  are  summed  in  M  successive  interpulse  intervals  to  produce  a 
detection  statistic  of  the  form 

M 

z“1/2^2rj,  (2.D 

j=l 

2 

where  rj  is  the  output  of  the  rectifier  in  the  j-th  interval.  When  suitably 
normalized,  Z  has  a  chi-squared  distribution  with  2M  degrees  of  freedom  under 
hypothesis  Hq  that  no  signal  is  present;  its  moment  generating  function  (m.g.f.) 
is 

hQ(u)  -  E(e'uZ|  HQ)  -  (1  +  u)'M.  (2.2) 

Under  hypothesis  Hj  that  a  target  is  present,  an  echo  is  received  in  the 
j-th  interval  with  energy  Ej  and  signal-to-noise  ratio  d^  *  2Ej/N,  1  <  j  <  M, 
where  N  is  the  unilateral  spectral  density  of  the  white-noise  background.  Then 
the  m.g.f.  of  Z  is 

hjCu)  «  E(e  1^)  =  (1  +  u)-*1  exp  (-  , 

,  rn  -  (2.3) 

S  -  l/2  D  -  V2  dk  ”  VN’ 

k-1 

where  Ep  is  the  total  received  energy.  The  probability  of  detection  is  given 
by  the  M-th  order  Q-function, 

9 

Qd  -  Pr  (Z  >  ZQ|  Hj)  -  QM(D,  m (2.4) 
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where  Zq  is  the  decision  level  with  which  the  statistic  Z  is  compared  [3,  pp. 
215-219].  The  false-alarm  probability 

Qo  '  qm(°*  ^  '  (2‘5) 

can  be  expressed  in  terms  of  the  incomplete  gamma  function.  A  saddlepoint 
approximation  was  applied  in  [4]  to  the  calculation  of  the  decision  level  Zq  and 
the  probability  of  detection  for  large  numbers  M  of  pulses  integrated;  typical 
results  are  tabulated  there  and  compared  with  the  exact  probabilities. 

When  the  radar  echoes  fade,  the  distribution  of  the  statistic  Z  must  be 
averaged  with  respect  to  the  distribution  of  signal  amplitudes.  It  is  simplest 
to  average  the  m.g.f.  in  (2.3)  and  invert  it  to  calculate  the  cumulative 
distribution  of  Z  and  thence  the  average  probability  of  detection.  Four  typical 
distributions  of  the  fading  amplitudes  were  enounced  by  Swerling  [5],  and  we  list 
them  here  with  the  resulting  m.g.f. 's  of  Z: 


Case  I:  D  *>  ADq,  where  A  has  a  Rayleigh  distribution, 

P(A)  =  (A/Aq2)  exp  (-A2/2Aq2),  A  >  0, 

yields  the  m.g.f. 

h^u)  -  (1  +  u)"(M_1)[l  +  (Bj  +  Du]*1 

-  V2  <D2)  =  Rj/N, 


(2.6) 


(2.7) 


(^indicating  an  average  with  respect  to  the  distribution  of  signal  amplitudes. 


Case  II:  dj  ■  AjdQ,  Vj,  where  the  Aj  are  statistically  independent  and  have 
Rayleigh  distributions  like  that  in  (2.6),  yields  the  m.g.f. 


7 


h2(u)  =  {1  +  (1  +  B2)u]  , 

B2  -  V2  <dj2>  «  (E.j>/N  - 


(2.8) 


The  distribution  of  the  statistic  Z  is  a  scaled  chi-squared  distribution  with  2M 
degrees  of  freedom. 


Case  III:  D  *  ADn,  where  A  has  the  density  function 


P(A)  -  (A3/2Aq2)  exp  (-A2/2A02),  A  >  0, 


yields  the  m.g.f. 


h3(u)  =  (1  +u)“(M"2)[l  +(1  +B3)u]-2, 
B3  -  I/4  <D2}  -  <Et)/2N. 


(2.9) 


(2.10) 


Case  IV:  dj  •»  AjdQ,  Vj,  where  the  Aj  are  statistically  independent  and  have  the 
distribution  of  Case  III,  yields  the  m.g.f. 


h4(u)  -  (1  +  u)M[l  +  [1  +  B4)u]'2M, 

**  V4  (d j2 )  -  (Ej  )/2N  -  B3/M. 


(2.11) 


The  distributions  of  the  sum  Z  arising  from  these  "Swerllng  cases"  are  given 
by  di  Franco  and  Rubin  [6,  ch.  11].  For  cases  III  and  IV  the  cumulative 
distributions  are  quite  complicated,  involving  hypergeometric  functions  of  orders 
depending  on  H.  For  case  IV,  for  instance,  calculating  the  detection  probability 
Qj  requires  summing  M  terms,  each  containing  an  incomplete  gamma-function  of 
order  M  +  k-1,  0  <  k  <  M. 

In  general,  if  -  . 
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G(u)  =  E(e"Su)  =  /  p(S)  e"uS  dS 


(2.12) 


Is  the  m.g.f.  of  the  positive-valued  random  variable  S,  representing,  the  total 
strength  of  the  fading  echoes,  the  m.g.f.  of  the  random  variable  Z  in  (2.1)  is, 
by  (2.3), 


h(u)  =  (1  +  u)"M  G  (y-^)  . 


(2.13) 


Its  singularities  will  always  lie  in  the  left  halfplane.  In  particular,  if  the 
signal  strength  has  a  gamma  distribution,  of  which  Sterling's  cases  I  -  IV  are 
special  instances. 


p(S)  =  [r(k)]_1  (k/S)k  Sk_1  e  ks/s, 


(2.14) 


where  S  is  the  mean  total  signal-to-noise  ratio,  then 


Su  •v-k 


C(»)  =  (1  +  f ) 


(2.15) 


and  the  m.g.f.  of  Z  is 


h(u)  =  (1  +  u)k^  (1  +  bu)  k,  b  =  1  +  S/k. 


(2.16) 


For  cases  I  -  IV,  k  -  1 ,  M,  2,  2M  respectively. 

The  method  of  integration  along  a  parabolic  contour,  as  in  (1.11),  has  been 
successfully  applied  to  evaluating  the  tail  probabilities  for  both  unfading 
signals  and  signals  having  random  strengths  described  by  distributions  p(S)  of 
the  form  in  (2.14).  A  bound  has  been  derived  for  the  truncation  error  incurred 
by  cutting  off  the  integration  in  (1.11)  at  a  finite  value  of  y.  In  most 
Instances  it  suffices  to  stop  the  numerical  integration  when  the  absolute  value 
of  the  Integrand  falls  below  a  fraction  e  of  the  accumulated  sum  times  the  step 
size  Ay,  whereupon  the  relative  truncation  error  will  be  less  than  e.  The  number 
of  steps  of  numerical  integration  required  to  attain  a  certain  accuracy  is  more 


or  less  independent  of  the  number  M  of  signals  processed  by  the  radar  receiver. 


(b)  Detection  in  Noise  of  Unknown  Level 

When  the  radar  is  being  jammed  by  a  transmitter  of  Gaussian  noise  spread 
over  a  frequency  band  much  wider  than  that  of  the  target  echoes,  the  receiver  is 

faced  with  the  problem  of  detecting  them  in  noise  of  unknown  and  random  spectral 

^  2 

density.  A  strategem  that  is  often  adopted  utilizes  K  samples  rj  of  the 
quadratically  rectified  output  of  a  narrowband  pass  filter  tuned  to  the  same 
frequency  as  that  of  the  echoes,  or  to  a  nearby  frequency;  these  samples  are 
taken  at  times  when  no  signals  are  expected  to  be  present,  and  it  is  presumed 
that  they  have  the  same  probability  distribution  as  the  noise  components  of  the 
terms  r^  in  (2.1).  The  strength  of  that  noise  is  then,  within  an  inessential 
constant  of  proportionality,  estimated  by  the  value  of  the  statistic 


r  E r;2- 

j=i 


(2.17) 


When  these  samples  r j  indeed  represent  noise  alone,  the  statistic  Z'  has  a 
scaled  chi-squared  distribution  with  2K  degrees  of  freedom  and  is  independent  of 
Z  in  (2.1).  The  decision  level  with  which  the  statistic  Z  is  compared  in  order 
to  decide  whether  any  radar  echo  is  present  is  made  proportional  to  Z',  or 
equivalently,  the  receiver  forms  the  statistic 

X  =  Z  -  SZ' ,  (2.18) 

and  it  decides  that  a  signal  is  present  whenever  X  >  0.  The  constant  p  is 
selected  to  achieve  a  pre-assigned  false-alarm  probability 


Q0(p)  =  Pr  (X  >  0|  HQ), 


(2.19) 


w* 


and  one  wishes  to  calculate  the  probability 

Qd< B;  S)  =  Pr  (X  >  0|  Hx)  .  -  (2.20) 

of  detecting  a  sequence  of  M  radar  echoes  of  total  signal-to-noise  ratio  S. 

If  the  signal  to  be  detected  has  a  fixed,  known  strength,  as  in  (2.3),  the 
statistic  X  has  the  m.g.f. 

n(u)  =  E(e  UX)  =  (1  +' u)  M(1  -  8u)  K  exp  (-  y )  (2.21) 

with  S  the  signal-to-noise  ratio  defined  in  (2.3).  The  probability  of  detection 
is  then  given  by 

Qd(6;  S)  =/  (1  +  u)-M(l  -  6u)_K  exp  (--p^-)  1^-,  (2.22) 

JC_ 

where  C_  is  a  contour  parallel  to  the  imaginary  u-axis  and  passing  between  the 
origin  and  the  essential  singularity  of  the  integrand  at  u  =  -1.  One  evaluates 
(2.22)  when  the  signal  strength  S  is  so  small  that  the  expected  value  of 
X  *  Z  -gZ'  is  negative;  here 


E(X)  =  S  +  M  -  8K. 


(2.23) 


When  E(X)  >  0,  on  the  other  hand,  one  evaluates 


i  -  V&;  s)  =  /(1+  U)’M(1  -  6u)"k  exP  (-  r^r)  Sr' 

along  a  contour  lying  in  0  <  Re  u  <  1/?. 

The  prescription,  "Choose  the  alternative  hypothesis  Hj  when  X  * 
Z  -gZ'  >  0,"  is  equivalent  to  the  F-test  of  the  a.nalysis  of  variance, 


(2.24) 


and  the  F- 


atatistic  is 
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F  -  KZ/MZ'. 


(2.25) 


The  false-alarm  probability  Qq( 3)  is  related  to  the  central  F-distribution,  and 
the  detection  probability  Q q(3;  S)  to  the  noncentral  F-distribution,  the  signal- 
to-noise  ratio  S  playing  the  role  of  the  noncentrality  parameter.  In  [7]  the 
computation  of  the  noncentral  F-distribution  by  integrating  (2.22)  or  (2. 24) 
numerically  along  a  straight  vertical  contour  was  treated,  and  a  bound  on  the 
error  incurred  by  truncating  the  range  of  numerical  integration  was  presented. 

It  was  shown  that  the  simple  rule  that  the  summation  be  stopped  when  the  absolute 
value  of  the  integrand  falls  below  a  fraction  e  of  the  accumulated  sum  tines  the 
step  size  Ay  suffices  to  bound  the  relative  error  in  the  computed  probability 
within  that  same  fraction  c. 

If  the  radar  echoes  are  fading  as  described  in  part  (a),  the  average 
detection  probability  is  obtained  by  averaging  the  integrals  in  (2.22)  and  (2.24) 
with  respect  to  S, 


•/ 


V6i  5)  -  /  <‘  +  -  *“>"*  G(rf^  75"’ 


i  -  Q. 


t(e;  sy  -  / 


(1  +  u)'M(l 


„  .-K  _/  u  i  du 
3u)  G(l  +  u  2  ni* 


(2.26) 


(2.27) 


where  5  is  the  average  total  signal-to-noise  ratio  and  G(u)  is  the  m.g.f.  of  the 
distribution  of  the  random  signal-to-noise  ratio  S,  as  in  (2.12).  Again  it  is 
possible  to  evaluate  these  probabilities  by  numerical  integration  along  a 
vertical  straight  line  or  other  suitable  contour.  Bounds  on  the  truncation  error 
when  the  m.g.f.  is  given  as  in  (2.15)  have  been  determined  and  will  be  presented 
In  a  paper  under  preparation.  p 
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t 


3.  The  Average  Power  of  a  Random  Process 


(3.1) 


z(t)  =  x(t)  +  s(t),  0  <  t  <  T,  (3.2) 

in  which  x(t)  is  a  stationary  Gaussian  random  process  of  mean  zero  and 
autocovariance  function 


We  seek  the  probability  distribution  of  the  average  power 


V  -  T_1  J  ( z(  t 


)]2  dt 


of  a  Gaussian  random  process 


4>(t  -  s)  *  Elx(t)  x(s)] , 


(3.3) 


and  s(t)  is  a  deterministic  signal.  This  problem  was  first  treated  by  Slepian 
[8],  who  assumed  s(t)  *  0;  references  to  other  work  are  given  in  the  proposal 
{9].  Our  approach  is  to  calculate  the  cumulative  distribution  q~(x)  =■  Pr  (V  <  x) 
and  its  complement  q+(x)  *  Pr  (V  >  x)  by  numerical  integration  of  the  contour 
Integrals 


q  (x) 


h(u)  e1^  du/2iti, 


q+(x)  a  -  J  u  *  h(u)  eUX  du/2iri, 


(3.4) 


(3.5) 


where  C_  and  C+  are  suitable  contours  obtained  by  deforming  the  straight  vertical 
contours  in  (1.5)  and  (1.6).  The  contours  pass  to  the  right  of  any  singularities 
of  the  moment-generating  function  (m.g.f.)  h(u),  which  in  this  problem  lie  on  the 
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negative  real  axis  of  the  u-plane.  The  contour  C_  lies  to  the  left  of  the 
origin,  C+  to  the  right. 

The  moment-generating  function  — 

_„v 

h(u)  =  E(e  )  (3.6) 

of  the  random  variable  V  was  shown  in  the  proposal  [9]  to  be  given  by 

h(u)  =  (D(2u/T)]_1/2  exp  [-  |^  J(2u/T)],  (3.7) 

where  D(u)  is  the  Fredholm  determinant 


D(u) 


n  ^  +  xju) 

j=i 


associated  with  the  homogeneous  integral  equation 


Af(t) 


4>( t  -  s)  f(s)  ds, 


(3.8) 


(3.9) 


whose  kernel  is  the  autocovariance  function  $(t  -  s)  and  whose  eigenvalues  A  are 
denoted  by  A^ ,  A^,  ...  .  The  function  J(u)  is  given  by 


where 


J(u) 


V 


T 

s(t)  f^(t)  dt 


(3.10) 


(3.11) 


are  the  coefficients  of  a  Fourier  series  for  the  signal  s(t)  in  terms  of  the 
eigenfunctions  fj(t)  of  (3.9).  As  indicated  in  the  proposal  [9], 


1A 


(3.12) 


c 

D(u)  “  exp  /  h(t,  t;  u;  t)  dt 
Jt\ 


in  terms  of  the  solution  h(r,  t;  u;  t)  of  the  integral  equation 


r 

h(r,  t;  u;  t)  +  u  /  $(r  -  v)  h(v,  t;  u; 


t)  dv  -  u$(r  -  t). 


(3.13) 


0  <  (r,  t)  <  x, 


J(u)  -  1/2  /  [s(t)  -  o(t;  u)]  dt 


r 

o(t;  u)  ■=  /  h(t,  v;  u;  t)  s(v) 


(3.14) 


(3.15) 


(a)  Application  to  RLC  Noise 

In  order  to  test  various  methods  of  calculating  the  m.g.f.  h(u)  needed  in 
(3.4)  and  (3.5)  and  of  carrying  out  the  subsequent  numerical  evaluation  of  those 
contour  integrals,  we  took  the  process  z(t)  *  x(t)  to  be  Gaussian  RLC  noise  of 


(3.16) 


spectral  density 

S(u>)  -  —2  2  2  2  2 

-  Uq  J  +4  p  a) 

and  unit  variance,  for  which  the  autocovariance  function  is 


♦(8)  -  (r^  -  i^)  [nj  exp  (— |  s  |  >  -  exp  (-t^M)],  (3.17) 


.  ,  2  2.1/2  ,  2  2  A/2 

hj  ■  U  +  (U  ~  i»Q  /  .  -  v  -  (p  -  ft>Q  ) 


(3.18) 


The  signal  s(t)  was  assumed  zero.  For  this  process  the  m.g.f.  h(u)  was  worked 


-i  3 

?  J 


r 


with  Q  the  usual  quality  factor  of  the  RLC  circuit,  Q  *»  u^/2ii. 

The  cumulative  distributions  for  the  average  power  V  of  this  RLC  noise  were 
computed  for  various  values  of  r  =  uT  and  Q  and  are  plotted  in  Figs.  1  and  2. 

The  contours  C+  and  C_  were,  as  discussed  in  Sec.  1,  taken  as  parabolas  passing 
through  saddlepoints  ug+  and  uq“,  respectively;  uq+  >  0,  and  uq~  <  0,  uq~  lying 
to  the  right  of  the  leftmost  singularity  of  the  m.g.f.  h(u).  The  integral  in 
(1.11)  was  evaluated  by  the  trapezoidal  rule.  The  size  of  the  segments  Ay  into 
which  the  integration  variable  y  is  divided  was  taken  initially  as 

Ay  *  [2*"(u0>]'1/2f  (3.20) 

and  it  was  successively  halved  until  the  value  of  the  integral  ceased  changing  in 
the  number  of  significant  figures  sought.  The  integration  was  stopped  when  the 
absolute  value  of  the  integrand  fell  below  10”^  times  Ay  times  the  accumulated 
sum. 

The  method  used  here  has  the  advantage  over  Slepian's  [8]  of  not  requiring 
prior  computation  of  the  eigenvalues  X  of  (2.9);  for  a  large  time-bandwidth 
product  r  ■  yT  the  number  of  significant  eigenvalues  X  is  of  the  order  of  r  and 
also  large,  and  Slepian's  method  entails  an  overlong  computation.  In  ours  the 
number  of  steps  of  numerical  integration  remains  of  roughly  the  same  order  of 


magnitude  over  the  entire  range  of  values  of  r  ■  pT. 

The  results  plotted  in  Figs.  1  and  2  show  that  the  distribution  of  V  is 
Insensitive  to  the  value  of  Q  unless  r  *  pT  is  of  the  order  of  1  or  less.  For 
r  **  4  and  r  ■  8  the  curves  for  Q  =  0  and  Q  =  »  lie  so  close  that  curves  for 
intermediate  values  of  Q  could  not  be  exhibited. 


(b)  Krein-Levinson  Method 

When,  as  we  assume  in  this  report,  the  noise  process  x(t)  is  stationary,  the 
functions  D(u)  and  J(u)  needed  for  the  m.g.f.  in  (3.7)  can  in  principle  be 
calculated  by  the  KreJn-Levinson  method  described  by  Kailath  et  al .  [10],  and  we 
have  begun  to  investigate  this  possibility. 

We  define 


Au(t,  r)  =  h(x,  r;  u;  t) 


in  terms  of  the  solution  of  (3.13).  It  obeys  the  integral  equation 


(3.21) 


r 

A>u(t,  r)  +  u  I  A^(x,  v)  4»( v  -  r)  dv  -  u<|>(  x  -  r) , 


which  is  obtained  by  setting  t  ■  x  in  (3.13)  and  using 


h(r,  t;  u;  x)  -  h(t,  r;  u;  x) , 


In  terms  of  this  function,  (3.12)  and  (3.15)  become 


0  <  r  <  x,  (3.22) 


D(u)  «  exp 


A- 


t)  dt, 


r 

o(t;  u)  ■  J  Ay(t,  v)  s(v)  dv. 


(3.23) 


(3.24) 


As  shown  in  [10],  this  function  Au(t,  r)  obeys  the  Sobolev  identity 


(It  +  If)  Au(t»  r)  “  _AuCt‘  0)  Au(t;  '  r)* 


(3.25) 


and  that  paper  suggests  solving  (3.22)  by  a  numerical  method  based  on  this 
identity  and  on  the  equation  (3.22)  with  r  ■  0, 


k  (t,  0)  +  u  I  A^(t,  v)  $(v)  dv  *=  u$(t). 
Jr\ 


(3.26) 


Taking  small  steps  of  size  A  in  the  variables  t  and  r,  one  writes  these 
approximately  as 


A  CE+Ta,  j+1  A)  =  A  (kA,  jA)  -  AA  (kA,  0)A  (kA,  k-jA) 
u  u  u  u 


and 


Au(kA,  0)  =  u$(kA)  -  uA  ^  A^kA,  jA)  $(jA). 

j-1 

Starting  with  the  initial  condition 


(3.27) 


(3.28) 


Au(0,  0)  -  u<K0)  “  u, 


(3.29) 


one  determines  AU(A,  A)  from  (3.27)  and  then  Au(A,  0)  from  (3.28).  Continuing, 
one  has  at  stage  k  the  values  of  Ay(kA,  jA),  0  <  j  <  k,  and  uses  them  in  (3.27) 
to  compute  Au(k+1A,  j+1 A),  0  <  j  <  k.  These  are  then  used  in  (3.28),  with  k 
replaced  by  k+1 ,  to  determine  A^(k+1 A,  0).  This  method  goes  by  the  name  of  the 
"KreJn-Levinson  algorithm”  [10].  The  integrals  in  (3.23),  (3.24),  and  (3.14), 
replaced  by  numerical  quadrature  formulas,  can  be  used  simultaneously  to  compute 
D(u)  and  J(u).  This  whole  computation  is  carried  out  for  each  point  u  on  the 
contour  of  Integration,  which  we  are  taking  to  be  the  parabola  figuring  in 
(1.11). 

This  algorithm  was  tried  for  RLC  noise,  and  it  was  found  inaccurate  even 
when  the  interval  (0,  T)  was  divided  into  one  hundred  steps,  A  “  0.01T.  The 
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T 


approximation  of  (3.25)  by  (3.27)  is  simply  too  crude.  Improved  accuracy,  yet 
with  fewer  subintervals  A,  resulted  from  replacing  the  integral  equation  (3.22) 
by  its  approximation  by  the  trapezoidal  rule. 


Au(kA,  iA)  +  uA  |  Ay(kA,  0)  4>(iA) 
k-1  L 

+  ^  Ay( kA,  jA) 

J=1 _ 

-  u  *(k-iA)  (3.30) 


♦(i-JA)  +  y 


A  (kA, 
u 


kA)  ^  k- 


-i  A)J 


which  constitutes  a  set  of  (k  +  1)  linear  simultaneous  equations  for 

A^(kA,  iA),  0  <  i  <  k.  These  were  solved  approximately  by  three  iterations  of 

the  conjugate-gradient  method  [11],  starting  with  an  approximate  solution 
obtained  from  (3.27)  and  (3.28). 

In  order  to  evaluate  the  contour  integral  (1.11)  by  numerical  integration 
along  a  parabolic  contour  through  the  saddlepoint  Uq,  it  is  necessary  first  to 
find  the  saddlepoint  by  solving  (1.8)  for  u  =  uQ  and  then  to  calculate  the 
curvature  k  from  (1.10),  and  this  entails  computing  the  first  three  derivatives 
of  the  phase  *(u)  of  (1.7)  with  respect  to  u  and  evaluating  them  at  the 
saddlepoint  uq,  which  is  real.  In  order  to  avoid  the  iteration  involved  in 
solving  (1.8)  for  u  **  uq,  we  picked  the  value  of  uq  at  the  start  and  used  (1.8) 

to  determine  the  corresponding  value  of  the  level  V  =  x  at  which  the  cumulative 

distribution  or  its  complement  is  evaluated. 

From  (3.7),  (1.7),  and  (3.23)  the  phase  of  the  integrand  is 


•(u)  *  D(2u/T)  +  ux  -  In  (±  u) 
»T 


■/ 


A^,(t,  t)  dt  +  ux  -  In  (±  u) ,  u'  *  2u/T. 


(3.31) 


The  derivatives  of  the  phase  $(u),  therefore,  involve  the  partial  derivatives  of 
A^(t,  r)  with  respect  to  the  real  parameter  u,  and  these  partial  derivatives  of 


AyCt,  r)  obey  integral  equations  obtained  by  differentiating  (3.22)  with  respect 
to  u.  The  Sobolev  identity  (3.25)  can  likewise  be  differentiated  with  respect  to 
u.  The  resulting  integral  equations  for  the  partial  derivatives  of  -A^t;  r),  u  ■ 
Uq,  were  also  solved  by  the  conjugate-gradient  method,  using  the  derivatives  of 
(3.27)  and  (3.28)  with  respect  to  u  to  provide  starting  values.  Because  the 
second  and  third  derivatives  are  not  needed  to  great  accuracy,  fewer  iterations 
of  the  conjugate-gradient  algorithm  suffice. 

This  scheme  succeeded  quite  satisfactorily,  but  further  study  and  refinement 
are  necessary.  The  solution  ^^(t,  r)  of  (3.22)  contains  a  term  of  the  form 
u(l  +  tu)  \  0  <  r  <  t,  which  dominates  when  pr  «  1 ,  p  the  damping  constant  in 
(3.16-.18).  When  the  imaginary  part  of  u  is  large,  that  is,  when  the  point  u  is 
far  out  on  the  parabolic  contour  of  integration,  this  term  results  in  a  rapid 
initial  gyration  of  the  solution  as  a  function  of  t,  and  this  cannot  be 
accurately  followed  by  the  numerical  algorithm  unless  very  small  step  sizes  A  are 
utilized.  The  ensuing  inaccuracy  is  not  serious,  for  the  value  of  the  integrand 
of  (1.11)  is  by  then  quite  small.  Nevertheless,  an  attempt  will  be  made  to 
subtract  this  term  from  A^( t,  r)  at  the  beginning  and  to  solve  only  for  the 
residual  portion  of  that  function  numerically. 

(c)  Chandrasekhar  Equations 

Alternative  methods  of  computing  the  m.g.f.  h(u)  are  available  when  the 
process  x(t)  can  be  considered  as  the  output  of  a  linear  system  driven  by  white 
noise.  Under  our  continuing  assumption  that  the  process  x(t)  is  stationary,  the 
system  can  be  taken  as  time  invariant  and  to  be  in  the  steady  state.  It  is  then 
necessary  to  know  an  n  x  n  dynamical  matrix  £,  an  n-element  column  vector  £,  and 
an  n-element  row  vector  £  such  that  the  process  can  be  represented  as 

x(t)  -  C  x(t)  (3.32) 

a# 
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with 


dx 

-  dT  “  —  +  £  w(t)»  (3*33) 

where  jc  is  an  n-element  state  vector  and  w(t)  is  white  Gaussian  noise  with 
spectral  density  R.  The  steady-state  variance  equation 

FA  +  AF+  +  R  GG+  =  0  (3.34) 

is  solved  for  the  steady-state  variance  matrix  and  the  known  value  of 

<f>(0)  =  CAC+  (3.35) 

establishes  the  required  value  of  R. 

As  indicated  in  the  proposal  [9],  the  Fredholm  determinant  D(u)  can  be 

calculated  from  the  solution  £  of  the  matrix  Riccati  equation 

~u 

d£u  +  +  + 

ir  “  ESu  +  -  S„£  ££u+  u  RS£  >  0  <  t  <  T,  (3.36) 

with  initial  condition 


£  (0)  «  uA(0); 

•MI  ~ 


(3.37) 


u  is  in  general  a  complex  variable  and  stands  for  a  point  on  the  contour  of 
integration.  Then 


D(u)  »  exp 


Furthermore,  the  function  o(t;  u)  in  (3.14)  is  given  by 

«Kt;  u)  »  ^(t),  . 


[/ 


C+  dt  . 


(3.38) 


(3.39) 
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where  fty(t)  is  the  solution  of  the  vector  differential  equation 

d§ 

dT”  Uu  +  C+[s(t)  -  o< t;  u)]  (3.40) 

in  terms  of  the  solution  of  (3.36). 

The  RLC  noise  treated  in  part  (a)  can  be  generated  by  a  second-order  system 
with  dynamical  matrix 


F  = 


(3.41) 


and  with 


G  « 


» 


C  -  (1  0). 


(3.42) 


The  matrix  is  then  ordinarily  2  x  2,  but  when  as  here  u  is  complex  and  one 
wishes  to  use  standard  routines  for  solving  systems  of  first-order  differential 
equations,  it  is  necessary  to  consider  the  real  and  imaginary  parts  of  the  four 
elements  of  £  (t)separately,  and  six  first-order  differential  equations  for  the 

Ml 

distinct  elements  of  this  complex  matrix  need  to  be  solved.  Simultaneously, 
(3.38)  is  integrated  after  conversion  to  a  differential  equation, 

~~  C£  (t)  C*  0  <  t  <  T,  E(0)  **  0,  (3.43) 

dt  — u  ~ 


D(u)  -  exp  E(T).  (3.44) 

Since  E(t)  is  also  complex,  this  adds  two  more  equations,  making  eight  in  all. 

In  general,  for  an  n-th  order  system,  the  number  of  differential  equations  to  be 
integrated  amounts  to  tr  +  n  +  2. 

By  using  the  DVERK  routine  in  the  IMSL  Library  of  computer  routines  we  were 
able  to  Integrate  the  Rlccati  equation  (3.36)  along  with  (3.43),  and  we  obtained 

r 

values  of  the  Fredholm  determinant  D(u)  in  agreement  with  those  calc elated  from 


22 


the  Chandrasekhar  method  requires  solving  An  +  2  first-order  differential 
equations,  which  amount  for  a  second-order  system  to  ten  instead  of  eight  as  with 
the  matrix  Riccati  equation.  For  a  third-order  system  both  methods  require 
solving  fourteen  equations,  and  for  n  >  3  the  Chandrasekhar  method  involves  fewer 
equations  than  the  Riccati  method. 

When  we  tried  the  Chandrasekhar  method  for  RLC  noise,  instability  arose  only 

when  the  imaginary  part  of  the  complex  parameter  u  was  large,  that  is,  at  points 

far  out  on  the  parabolic  contour  of  integration  in  (1.11).  This  appears  to  be 

related  to  the  term  u(l  +  xu)“*  that  causes  an  initial  sharp  gyration  in  the 

* 

function  A  (x,  r)  treated  in  part  (b).  Indeed,  A  (x,  x)  *>  C  K  (x),  0  <  x  <  T. 

u  U  ~  "il 
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We  shall  have  to  Investigate  whether  it  is  possible  to  subtract  out  what 
corresponds  to  that  term  in  order  to  make  the  Chandrasekhar  method  also  efficient 
for  all  values  of  u  at  which  the  function  D(u)  is  needed.  — 

Since  (3.40)  depends  on  only  through  =  £^C+,  the  Chandrasekhar  method 
will  also  suffice  for  calculating  the  function  J(u)  in  (2.7)  and  hence  can  be 
applied  even  when  the  process  z(t)  contains  a  deterministic  signal  s(t),  but  this 
has  not  yet  been  studied. 

The  derivatives  of  the  phase  $(u)  of  the  integrands  of  (3.4)  and  (3.5), 
defined  as  in  (1.7),  can  also  be  calculated  by  the  Chandrasekhar  method.  The 
necessary  equations  are  obtained  by  differentiating  (3.45),  (3.46),  and  (3.48) 
with  respect  to  u.  These  derivatives  are  needed  only  for  real  values  of  u.  For 
RLC  noise  we  calculated  the  first  three  derivatives  of  4>(u)  at  the  saddlepoint  by 
integrating  the  fifteen  simultaneous  first-order  differential  equations  that 
arise  by  this  process,  using  them  as  described  at  the  end  of  part  (b)  to 
determine  the  level  V  =  x  by  (1.8)  and  the  curvature  k  of  the  path  of  integration 
by  (1.10).  The  results  agreed  well  with  those  calculated  by  the  exact  m.g.f.  in 
(3.19).  A  detailed  comparison  of  the  methods  of  parts  (b)  and  (c)  remains  to  be 
carried  out.  The  former  requires  only  the  autocovariance  function  4>(t  -  r)  of 
the  process  x(t);  for  the  latter  a  state-space  model  must  be  constructed,  but 
standard  routines  for  integrating  systems  of  first-order  differential  equations 
can  then  be  applied. 
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4.  Photoelectron  Counting  Probabilities 


(a)  Introduction 

Let  x  be  a  positive  integer-valued  random  variable  with  a  probability  dis¬ 
tribution 

Pr  (x  =  k)  =  p^,  k  =  0,  1,  2,  ..., 


cumulative  distribution  (c.d.) 

n-1 

On"  =  Pr  <x  <  R)  V 

k=0 

and  complementary  cumulative  distribution  (c.c.d.) 

00 

Qn+=  Pr  (x  >  n)  =  1  -  =  J^Pk. 

k=n 


(4.1) 


(4.2) 


In  the  applications  we  usually  have  in  mind,  x  is  the  number  of  electrons  counted 
during  a  specific  interval  (0,  T)  at  the  output  of  some  photoelectric  device. 

The  c.d.  and  the  c.c.d.  can  be  computed  from  the  probability  generating  function 
(p.g.f .) 


h(z)  =  E(z*)  =  ^  ^  pkzk 
k=0 


by  the  contour  integrals 


(4.3) 


(4.4) 


(4.5) 
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in  which  C_  and  C+  are  circles  centered  at  the  origin  z  *  0;  C+  also  encloses  the 
point  z=l,  C_  does  not;  and  the  contours  enclose  no  singularities  of  the  p.g.f. 
h(z),  all  of  which  lie  outside  the  unit  circle. 

As  shown  in  [13],  the  numerical  evaluation  of  these  integrals  is  most  effi¬ 
cient  when  the  contours  pass  through  saddlepoints  of  the  integrand,  that  is, 
through  points  Zq  that  are  solutions  of 


where 


*'(  z) 


h' ( z)  n  1 

h(  z)  z  z  -  1 


0, 


(4.6) 


T(z)  =  In  h( z)  -  n  In  z  -  ln[±(z  -  1)] 


(4.7) 


is  the  "phase"  of  the  integrands,  ±(z  -  1)  standing  for  1  -  z  when  (4. A)  is  inte¬ 
grated  and  for  z  -  1  when  (4.5)  is  integrated.  One  such  saddlepoint,  zq  =  zq, 
lies  in  0  <  Zq  <  1,  the  other  lies  on  the  real  z-axis  between  z  =  1  and  the  left¬ 
most  singularity  of  h(z).  The  saddlepoint  can  usually  be  most  expeditiously 
located  by  solving  (4.6)  by  Newton's  method  or  the  secant  method.  Alternatively, 
the  equivalent  of  (1.9)  can  be  used. 

The  circular  contours  in  (4.4)  and  (4.5)  often  suffer  the  disadvantage  that 
the  integrand  is  oscillatory  along  them,  so  that  too  many  steps  may  be  needed  in 
order  to  obtain  an  accurate  value  of  or  Q“  by  numerical  integration.  The 
value  of  the  integrand  decreases  most  rapidly  from  its  value  at  the  saddlepoint 
when  the  contour  is  taken  as  the  path  of  steepest  descent  through  the  saddle¬ 
point.  Along  that  path  Im  f(z)  =  0.  Determining  the  path  of  steepest  descent, 
however,  would  much  protract  the  evaluation  of  the  contour  integral.  In  many 
cases  it  has  been  found  that  the  path  of  steepest  descent,  which  lies  symmetri¬ 
cally  about  the  real  axis,  has  nearly  a  parabolic  form,  and  it  is  then  expedi¬ 
tious  to  use  instead  its  osculatory  parabola,  that  is,  a  parabola  passing  through 
the  saddlepoint  Zq  and  having  the  same  curvature 
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k  =  y  r'(z0)/r(z0)  (a. 8) 

at  the  point  z  ■  Zq.  (Primes  indicate  differentiation.)  Along  this  parabola 

2  =  zQ  +  1/2  <y^  +  iy,  (A. 9) 

and  the  integral  to  be  evaluated  numerically  is 

%*  "  "  7  116  [iff-  n  “  '  lK5,)]  d?-  <4-10> 

The  initial  step  size  is  taken  as 

Ay  -  ir(*0>r1/2  (A. ll) 

and  successively  halved  until  the  values  of  the  probability  stabilize  to  the 
number  of  significant  figures  desired.  The  trapezoidal  rule  is  used  for  reasons 
described  by  Rice  [2]. 


(b)  Single-Stage  Photomultiplication 

The  contour-integration  method  has  been  applied  to  finding  the  distribution 
of  the  number  of  electrons  at  the  output  of  a  photomultiplier  with  a  single  stage 
of  multiplication  [14].  Light  striking  a  photoelectrically  emissive  surface 
ejects  a  number  k  of  photoelectrons  during  an  interval  (0,  T)  with 
probability  11^.  The  p.g.f.  of  this  distribution  is 


f(z) 


(A. 12) 


Each  such  "primary"  electron  is  accelerated  by  an  applied  voltage,  impinges  on  a 
second  metallic  surface,  and  ejects  from  it  a  random  number  of  secondary 
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where 


(4.18) 


f 


k 


f(z) 


-G’ 

z=e 


and  where  S(n,  k)  are  modified  Stirling  numbers  of  the  second  kind  and  obey  the 
recurrent  relations 


S(1 ,  1)  =  1,  S(k,  n)  =  0,  k  >  n, 

S(n  +  1,  k)  =  [S(n,  k  -  1)  +  kS(n,  k)]/(n  +  1). 


(4.19) 


These  formulas  were  used  to  compute  exact  values  of  the  probabilities  and  Qn 
for  comparison  with  those  computed  by  our  approximation  methods.  They  involve 
lengthy  computations  when  the  number  n  is  large,  and  round-off  error  introduces 

when  is  computed  as  in 

(4.1). 

Three  types  of  primary  distribution  {n^}  were  treated:  (i)  that  arising 
when  the  incident  light  has  a  Lorentz  spectral  density,  and  for  which  the  proba¬ 
bilities  nk  can  be  calculated  by  a  method  given  by  Bddard  [15];  (ii)  the  nega¬ 
tive  binomial  distribution 


large  relative  error  into  the  c.c.d.  =  1  - 


\  -  -  V«V  +  » •  <4-20> 

with  Np  the  mean  number  of  primary  electrons,  and  M  the  number  of  temporal  modes 
or  "degrees  of  freedom" ,  given  roughly  by  the  product  of  the  bandwidth  W  of  the 
incident  light  and  the  duration  T  of  the  counting  interval;  and  (iii)  the  Poisson 
distribution 


1^  -  Npk  exp  (-Np)/k! ,  (4.21) 

for  which  the  output  distribution  possesses  the  Neyman  Type-A  form. 

For  primary  distributions  of  type  (i)  and,  for  M  an  integer,  of  type  (ii), 
Che  p.g.f.  f(z)  possesses  poles  z^  lying  on  the  positive  real  z-axis  to  the  right 
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A  . 


\ 


'•dr- 


of  z  *  1,  whereupon  the  output  p.g.f.  h(z)  possesses  poles  at  the  points 


4r)  =  1  +  G_1(ln  zk  +  2rirt),  (4.22) 

for  all  positive  and  negative  integral  values  of  r,  including  r  =  0.  The  contour 
integral  (3.5)  can  then  be  evaluated  as  a  residue  series  with  terms  corresponding 
to  each  point  of  this  doubly  infinite  lattice  of  poles.  For  the  negative-binomi¬ 
al  primary  distribution,  the  array  shrinks  to  a  vertical  column  of  multiple  poles 
of  order  M  at  the  points 

tr  =  1  +  G  * (-In  v  +  2rTri).  (4.23) 

When  the  number  n  of  output  electrons  is  of  the  order  of  or  greater  than  the 
expected  value 


E(n)  *  N  G,  (4.24) 

P 

only  a  very  few  of  the  poles  above  and  belcw  the  real  axis  contribute  signifi¬ 
cantly  to  the  residue  series.  The  number  of  columns  of  the  lattice  that  need  to 
be  included  under  case  (i)  (incident  light  with  a  Lorentz  spectral  density)  is  of 
the  order  of  the  time-bandwidth  product  WT.  Formulas  for  computing  the  residues 
and  numerical  examples  are  given  in  [14],  Unless  WT  is  very  large,  this  residue 
method  is  most  efficient  for  computing  the  c.c.d.  Qn+  with  n  >  E(n) . 

Consideration  was  given  to  using  the  contour  integrals  in  (4.4)  and  (4.5) 

for  computing  the  c.d.  and  the  c.c.d.  of  the  number  of  output  electrons.  The 

/ 

path  of  steepest  descent  is  now  too  complicated  to  be  utilized  or  even  approxima¬ 
ted,  for  it  consists  of  a  number  of  hairpin-like  curves  opening  to  the  right  and 
passing  around  each  of  the  poles  of  h(z).  We  therefore  used  instead  a  straight 
vertical  path  of  integration  passing  through  Che  saddlepoint  of  the  integrand 
lying  on  the  real  z-axis;  the  one  lying  in  0  <  t  <  1  was  used  for  ,  and  the 
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one  lying  in  1  <  z  <  was  used  for  Qn+»  these  having  been  determined  by  solv¬ 

ing  (4.6)  by  Newton's  method.  A  bound  on  the  error  incurred  by  truncating  the 
numerical  integration  was  worked  out.  Values  of  and  Q^+  over  a  broad  range 
of  values  of  n  could  be  accurately  computed  in  this  way  without  requiring  more 
than  one  hundred  or  so  steps  of  numerical  integration.  Details  are  given  in  Sec. 
4  of  [14].  This  method  does  not  require  that  the  number  K  in  (4.20)  be  an  inte¬ 
ger.  The  contour-integration  method  can  be  used  for  all  three  types  of  primary 
distribution. 

On  each  branch  of  the  path  of  steepest  descent  of  the  integrand  in  (4.4)  and 
(4.5)  there  is  a  point  zk  at  which  the  magnitude  of  the  integrand  is  maximum; 
this  is  a  saddlepoint  of  the  integrand.  The  contribution  of  each  branch  to  an 
integration  of  (4.4)  or  (4.5)  along  its  path  of  steepest  descent  can  be  crudely 
approximated  by  assuming  that  in  the  neighborhood  of  each  saddlepoint  zk  the 
integrand  has  a  Gaussian  form,  and  this  leads  to  the  formula 

OO 

Q  -  y;  [2itr(zk)]"1/2  exp  [y(zk)l,  (4.25) 

k=-«° 

where  T(zk)  is  the  phase  of  the  integrand  as  in  (4.7)  and  H"(zk)  is  its  second 
derivative  evaluated  at  the  k-th  saddlepoint  [4],  For  the  most  part  two  or  three 
saddlepoints  on  each  side  of  the  principal  one  Zq  on  the  real  z-axis  are  all  that 
it  is  worth  taking  into  account;  often  the  one  at  Zq  alone  suffices.  Improved 
accuracy  was  obtained  by  evaluating  the  contribution  of  the  principal  branch  of 
the  path  of  steepest  descent,  which  crosses  the  real  axis  at  zq,  by  means  of  the 
uniform  asymptotic  expansion  described  in  [16]  and  [17].  The  details  of  this 
method  are  to  be  found  in  Sec.  5  of  [14]. 


(c)  Multistage  Photomultiplier 


A  photomultiplier  usually  has  not  one,  but  several  stages  of  electron  multi 

plication.  We  number  the  stages,  each  associated  with  an  electrode  .called  a 

"dynode",  from  the  last.  A  single  electron  striking  the  k-th  dynode  ejects  n 

(k) 

secondary  electrons  from  it  with  probability  pR  ;  the  collection  of  these  proba 
bilitles  constitutes  the  k-th  "single-stage"  distribution  and  has  a  p.g.f. 


gk(z)  *"•  (4-26) 

n=0 

Then  if  there  are  N  stages  in  all,  the  probability  pn  that  n  electrons  emerge  at 
the  output  when  a  single  primary  electron  strikes  the  N-th  dynode  has  a  p.g.f. 

00 

Vz)  “XjPm2”*  (4,27) 

m*=0 

where 


Hl<z>  *  gj(z), 

\<z)  **  gk(Hk_j  (  z) ) »  k  =  2,  3,  ....  N. 


(4.28) 


If  the  distribution  of  primary  electrons  again  has  the  p.g.f.  f(z)  defined  in 
(4.12),  the  total  number  of  output  electrons  has  the  p.g.f.  f(HN<z)). 

In  this  study  we  have  concentrated  on  calculating  the  cumulative  output 
distributions,  defined  as  in  (4.1)  and  (4.2),  for  a  single  primary  electron,  so 
that  in  (4.3)-(4.5) 

h(z>  -  yo. 

We  assumed  that  all  the  single-stage  distributions  have  the  same  negative-binomi 
al,  or  PolyA  form,  with  common  p.g.f. 
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gk(z)  =  g(z)  “  [1  -  bG(z  -  1)]  ^b,  b  >  0,  (4.29) 

which  corresponds  to  the  distribution  in  (4.20)  with 

M  =  1/b,  v  =  bG/(l  +  bG). 

Again  G  is  the  gain  per  stage;  the  overall  gain  of  the  photomultiplier  is 

GQ  =  GN,  (4.30) 

where  N  is  the  number  of  stages.  The  Poisson  single-stage  distribution  represen¬ 
ted  by  (4.15)  arises  in  the  limit  b  ♦  0. 

The  singularities  of  the  overall  p.g.f.  H^(z)  now  lie  on  or  close  to  the 
portion  of  the  real  z-axis  lying  to  the  right  of  the  point  z  =  1 ,  and  the  path  of 
steepest  descent  has  roughly  the  form  of  a  parabola  lying  symmetrically  about  the 
real  z-axis  and  opening  to  the  right.  This  path  can  be  replaced  by  its  oscula- 
tory  parabola,  whose  curvature  k  at  the  saddlepoint  is  determined  by  (4.8).  The 
first  three  derivatives  of  the  phase  ¥(z)  of  the  integrand,  defined  as  in  (4.7), 

are  computed  by  an  N-fold  set  of  three  recurrent  relations.  The  cumulative  prob- 
“f* 

abilities  are  then  computed  by  numerical  integration  of  (4.10).  The  results 
were  compared  with  values  computed  by  the  recurrent  relations  given  by  Prescott 
[18]  for  N  =  3  and  values  of  n  up  to  100,  and  the  agreement  was  excellent.  The 
recurrent  relations  require  of  the  order  of  n  operations  to  compute  pn  or  , 
and  this  number  is  independent  of  the  number  of  significant  figures  sought.  Thus 
for  n  =  3000  of  the  order  of  nine  million  operations  would  be  required,  in  con¬ 
trast  to  the  ten  or  twenty  steps  of  numerical  integration  of  (4.10)  needed  for 
adequate  precision. 

When  the  cumulative  distributions  so  calculated  are  plotted  on  a  semi- 
logarithmic  grid  for  various  numbers  N  of  stages  of  multiplication,  the  overall 
gain  Gq  ■  GN  being  kept  fixed,  one  finds  the  points  lying  closer  and  closer  to 
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straight  lines  as  N  increases.  The  results  are  fitted  quite  closely  by  a 
geometric  or  "exponential"  distribution, 


1  -  (1  -  v)Gq, 

' 

a 

'c 

> 

CM 

CN 

> 

1 

V-/ 

o 

n  >  0, 

(A. 31) 

Gq(1  -  v)vn  1 , 

n  >  0, 

in  which  the  parameter  v  is  given  by 


v  *  a/(a  +1), 


*2<G0  -  l) 
2G(G  -  1)  * 


(A. 32) 


where  iCj  is  the  second  factorial  moment  of  the  single-stage  distribution  and  G  is 
the  gain  per  stage, 


G-g'(l),  K2-g”(l),  (A. 33) 

primes  denoting  differentiation.  This  value  of  v  is  chosen  so  that  the  distribu¬ 
tion  in  (A. 31)  has  the  same  variance  as  the  exact  distribution, 

°^G0^G0  ^  2  2 
Var  n  *  — g(g~-~T) — ’  ff  =  +  G  “  G  *  (A.3A) 

A  number  of  observers  have  reported  approximately  exponential  distributions 
for  the  number  of  output  electrons  in  photomultipliers,  but  an  adequate  explana¬ 
tion  of  this  phenomenon  seems  to  have  been  lacking  [18].  In  [19]  the  distribu¬ 
tion  in  (A. 31)  was  shown  to  approach  the  true  output  distribution  under  the  con¬ 
ditions  stated:  the  overall  gain  Gq  ■  GN  is  fixed,  the  single-stage  distribu¬ 
tions  are  identical,  and  the  number  N  of  stages  Increases,  the  gain  G  per  stage 
decreasing  simultaneously  toward  1 . 

This  behavior  is  the  counterpart  for  multiplicative  processes  of  the 

* 

central-limit  theorem  for  sums  of  Independent,  identically  distributed  random 
variables.  If  N  random  variables  have  mean  zero,  and  one  scales  their  sum  by 
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r 


-1  /2 

multiplying  by  N  ,  the  distribution  of  the  scaled  sura  approaches  the  Gaussian 
form  as  N  increases.  Furthermore,  if  the  random  variables  already  have  a 
Gaussian  distribution,  their  sum  has  exactly  a  Gaussian  distribution.  For  multi¬ 
plicative  processes  the  geometric  distribution  plays  the  same  role  as  the 
Gaussian  does  for  sums.  If  the  single-stage  distributions  are  geometric,  w4th 
p.g.f .'s 


gku> 


G.  (z  -  1) 

1  + - - - 

1  -  ak(z  -  1)* 


(4.35) 


the  output  distribution  is  also  geometric  as  in  (4.31),  with  p.g.f. 


Vz) 


G0(z  -  1) 

+  1  -  a(z  -  1)’ 


(4.36) 


in  which 


N 


k»l 


and  the  coefficients  afc.  The  operations 
in  (4.28)  then  constitute  a  sequence  of  homographic  transformations  of  the  com¬ 
plex  plane  that  leave  the  point  z  *=  1  invariant ,  and  such  transformations  form  a 
group.  One  might  expect,  therefore,  that  even  if  the  single-stage  distributions 
are  not  of  the  geometric  type,  the  output  distribution  would  approach  the  geomet¬ 
ric  form  in  (4.31)  as  the  number  of  stages  increases,  the  overall  gain  Gq  remain¬ 
ing  fixed,  and  this  behavior  was  demonstrated  in  [19].  Unlike  the  limiting  be¬ 
havior  predicted  by  the  central-limit  theorem  for  sums,  however,  the  variance  of 
the  output  distribution  of  electrons  continues  to  increase  with  an  increasing 
number  N  of  stages,  even  though  its  mean  Gq  -  GN  remains  fixed. 

Although  the  limiting  behavior  of  the  output  distribution  as  N  increases  has 
been  extensively  treated  in  the  literature  on  branching  processes,  the  emphasis 


and  a  depends  on  the  individual  gains  G^ 
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has  been  on  the  behavior  of  the  random  variable  W  =*  G  ^n  where  n  is  the  total 
number  of  output  electrons,  and  the  gain  G  per  stage  is  kept  fixed  in  the  passage 
to  the  limit  N  +  The  type  of  limiting  behavior  described  in  [19]- seems  not  to 
have  been  considered. 
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distribution  q  (x)  of  average  power  of 
indexed  with  the  value  of  Q. 


